smooth.max<-function(x,y,h=0.1) {
  (x*exp(h*x)+y*exp(h*y))/(exp(h*x)+exp(h*y))
}

u.pay<-function(p,r,inv,par.setting=NULL,par.behavior=NULL) {

  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  #old.payment.outcome<-par.setting[7]
  old.payment.outcome<-0 # set norm to not pay
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  period<-par.setting[11]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)

  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
  par.pay.chat<-0
  
  
#  (100-r)*p-inv.cost*inv-p*penalty+p*pay.bias+par.social.norm*(p*old.payment.outcome+(1-p)*(1-old.payment.outcome)) # both side social norm model
#  (100-r)*p-inv.cost*inv-p*penalty+p*pay.bias+(par.social.norm)*(p*old.payment.outcome+(1-p)*(1-old.payment.outcome)) # both side social norm model
#  (100-r)*p-inv.cost*inv-p*penalty+(par.social.norm+par.pay.chat*chat.payment)*(p*old.payment.outcome+(1-p)*(1-old.payment.outcome)) # both side social norm model

  
#  (100-r)*p-inv.cost*inv-p*penalty+(par.social.norm)*(p*old.payment.outcome+(1-p)*(1-old.payment.outcome))+par.pay.chat*chat.payment*(1-p) # both side social norm model
#  (100-r)*p-inv.cost*inv-p*penalty+(par.pay.social.norm+period*par.pay.social.learn)*((1-p)*(1-old.payment.outcome)) # both side social norm model
#  (100-r+par.pay.fairness)*p-inv.cost*inv-p*penalty+(par.pay.social.norm)*((1-p)*(1-old.payment.outcome)) # both side social norm model
   (100-r-par.pay.fairness*smooth.max(r-(100-r-inv.cost*inv),0))*p-inv.cost*inv-p*penalty+(par.pay.social.norm)*((1-p)) # both side social norm model
#  (100-r-par.pay.social.norm)*p-inv.cost*inv-p*penalty # both side social norm model
  
  
    #  (100-r)*p-inv.cost*inv-p*penalty+(par.social.norm+par.pay.chat*chat.payment)*((1-p)*(1-old.payment.outcome)) 
  
    #  (100-r)*p-inv.cost*inv-p*penalty+p*pay.bias
}

prob.pay<-function(r,inv,par.setting=NULL,par.behavior=NULL) {
  
  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
  up<-u.pay(1,r,inv,par.setting=par.setting,par.behavior=par.behavior)
  unp<-u.pay(0,r,inv,par.setting=par.setting,par.behavior=par.behavior)
  
#  sc.p<-gamma.pay*up+(par.social.norm)*old.payment.outcome # both side social norm model
#  sc.np<-gamma.pay*unp+(par.social.norm)*(1-old.payment.outcome) # both side social norm model

#  1/(1+exp(-(sc.p-sc.np)))  
  1/(1+exp(-gamma.pay*(up-unp)))  
}

u.r<-function(r,inv,par.setting=NULL,par.behavior=NULL) {

  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
#  r*(prob.pay(r,inv,par.setting=par.setting,par.behavior=par.behavior)+par.r.past.pay*old.payment.outcome)
  r*(prob.pay(r,inv,par.setting=par.setting,par.behavior=par.behavior))
}

qre.r<-function(inv,par.setting=NULL,par.behavior=NULL) {

  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
#  if(is.na(e0)) browser()
  r.list<-c(1:e0)
  
  ur<-u.r(r.list,inv,par.setting=par.setting,par.behavior=par.behavior)
  
  sc<-exp(gamma.r*(ur-max(ur)))
  pp<-sc/sum(sc)
  
#  browser()
  
  cbind(r.list,pp)
}

prob.success<-function(inv,par.setting=NULL,par.behavior=NULL) {
  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
  
  p0-p1*inv
}


u.att<-function(att,inv.list,par.setting=NULL,par.behavior=NULL) {
  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
#  result<-ea
  result<-ea
  if(att>0) {
    inv<-inv.list[att]
    
#    ps<-prob.success(inv,par.setting=par.setting,par.behavior = par.behavior)+prev.att*old.attack.outcome
    ps<-prob.success(inv,par.setting=par.setting,par.behavior = par.behavior)
    qr<-qre.r(inv,par.setting=par.setting,par.behavior=par.behavior)
    pp<-prob.pay(qr[,1],inv,par.setting=par.setting,par.behavior=par.behavior)
    
    result<-ps*sum(qr[,1]*qr[,2]*pp)
  }
  result
}

qre.att<-function(inv.list,par.setting=NULL,par.behavior=NULL) {

  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  inv.social.norm<-par.setting[10]
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.inv.chat<-0
  
  att.list<-c(0,1,2)
  ua<-NULL
  for(att in att.list) {
    ua<-c(ua,u.att(att,inv.list,par.setting = par.setting,par.behavior = par.behavior))
  }
  
  # flat QRE

  sc<-exp(gamma.a*(ua-max(ua)))
  pp<-sc/sum(sc)
  
  # nested QRE
  
  # p.d1<-1/(1+exp(-gamma.d*(ua[2]-ua[3])))
  # p.d2<-1/(1+exp(-gamma.d*(ua[3]-ua[2])))
  # 
  # u.att<-p.d1*ua[2]+p.d1*ua[3]
  # p.att<-1/(1+exp(-gamma.a*(u.att-ua[1])))
  # 
  # pp<-c(1-p.att,p.att*p.d1,p.att*p.d2)
  
#  if(pp[1]>0.99)  browser()
  
  cbind(att.list,pp)
}

u.inv<-function(inv,inv.other,inv.lag,par.setting = NULL,par.behavior = NULL) {
  
  # inv.lag = (me.lag, you.lag)

  e0<-par.setting[1] 
  ea<-par.setting[2]
  inv.cost<-par.setting[3]
  p0<-par.setting[4]
  p1<-par.setting[5]
  penalty<-par.setting[6]
  old.payment.outcome<-par.setting[7]
  chat.invest<-par.setting[8]
  chat.payment<-par.setting[9]
  # inv.social.norm<-par.setting[10]
  
  Period<-par.setting[11]
  
  inv.social.norm<-1
  
  # from llk.i   par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom)
  
  gamma.pay<-par.behavior[1]
  gamma.r<-par.behavior[2]
  gamma.a<-par.behavior[3]
  gamma.i<-par.behavior[4]
  par.me.lag<-par.behavior[5]
  par.pay.social.norm<-par.behavior[6]
  par.pay.fairness<-par.behavior[7]
  par.inv.social.norm<-par.behavior[8]
  par.you.lag<-par.behavior[9]
  
  par.r<-0
  par.inv.chat<-0
  
  
#  inv.cost<-inv.cost*(1+ps.bias)-(par.inv.social.norm+par.inv.social.learn*Period+par.inv.chat*chat.invest)*(inv*inv.social.norm)
#  inv.cost<-inv.cost*(1+ps.bias)-(par.inv.social.norm+par.inv.social.learn*Period)*(inv*inv.social.norm)
#  inv.cost<-inv.cost-(par.inv.social.norm)*(inv*inv.social.norm)
  
#  inv.cost<-inv.cost-par.me.lag*inv.lag[1]-par.you.lag*inv.lag[2]
  inv.cost<-inv.cost-par.me.lag*inv.lag[1]-par.you.lag*inv.lag[2]-(par.inv.social.norm)
  
  inv.list<-c(inv,inv.other)
  qra<-qre.att(inv.list,par.setting=par.setting,par.behavior=par.behavior)

  # bias of attack prob
    
  #qra[1,2]<-qra[1,2]*(1+pa.bias)
  #qra[,2]<-qra[,2]/sum(qra[,2])
  
#  browser()
  
  ps<-prob.success(inv,par.setting=par.setting,par.behavior = par.behavior)
#  ps<-p0*(1+ps.bias)-p1*inv
#  u.noloss<-(e0-inv*inv.cost)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps))+par.inv.ransom*(inv*inv.social.norm+(1-inv)*(1-inv.social.norm)) # change par ransom to inv social norm
#  u.noloss<-(e0-inv*inv.cost)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps))+(par.inv.social.norm)*(inv*inv.social.norm+(1-inv)*(1-inv.social.norm)) # change par ransom to inv social norm
#  u.noloss<-(e0-inv*inv.cost)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps))+(par.inv.social.norm)*(inv*inv.social.norm+(1-inv)*(1-inv.social.norm)) # change par ransom to inv social norm
  
#  u.noloss<-(e0-inv*inv.cost)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps))+
#    (par.inv.social.norm)*(inv*inv.social.norm+(1-inv)*(1-inv.social.norm))+par.inv.chat*chat.invest*inv # change par ransom to inv social norm

u.noloss<-u.risk(e0-inv*inv.cost,par.r)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps)) # change par ransom to inv social norm

#u.noloss<-(e0-inv*inv.cost)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps))+
#  (par.inv.social.norm)*(1+par.inv.chat*chat.invest)*(inv*inv.social.norm) # change par ransom to inv social norm

  
#    u.noloss<-(e0-inv*inv.cost)*(qra[1,2]+qra[3,2]+qra[2,2]*(1-ps))+
#    (par.inv.social.norm)*(inv*inv.social.norm+(1-inv)*(1-inv.social.norm))+
#    par.inv.chat*chat.invest*inv
  
  qr<-qre.r(inv,par.setting=par.setting,par.behavior=par.behavior)
  pp<-prob.pay(qr[,1],inv,par.setting=par.setting,par.behavior=par.behavior)
  
#  u.loss<-qra[2,2]*ps*sum(qr[,2]*(pp*(e0-qr[,1]-inv*inv.cost-penalty)+(1-pp)*(-inv*inv.cost)))  
  
#  rr<-(1-par.inv.ransom)*qr[,1]+par.inv.ransom*old.ransom
  rr<-qr[,1] # no anchoring
  u.loss<-qra[2,2]*ps*sum(qr[,2]*(pp*u.risk(e0-rr-inv*inv.cost-penalty,par.r)+(1-pp)*u.risk(-inv*inv.cost,par.r)))  
  
#  if(abs(u.loss+u.noloss)>100) browser()
#  if(gamma.i==0) browser()
 
#  if(is.na(u.loss+u.noloss)) browser() 
  #u.noloss+u.loss+par.inv.initia*(inv*inv.lag+(1-inv)*(1-inv.lag))
  
  result<-u.noloss+u.loss
  ifelse(result<(-500),500,result)
}

#u.risk<-function(x,r) {
#  ifelse(r==0,x,(1-exp(-r*x))/r)
#}

#u.risk<-function(x,r) {
#  ifelse(x>0,x,(1+r)*x)
#}

u.risk<-function(x,r) x

llk.i<-function(index,par,par.restrict=NULL,par.restrict.value=NULL,data=NULL) {
  if(!is.null(par.restrict)) {
    final.par<-par.restrict.value
    counter<-1
    for(i in 1:length(par.restrict)) if(par.restrict[i]==1) {
      final.par[i]<-par[counter] 
      counter<-counter+1
    }
  } else final.par<-par
  
#  if(index==1) cat(final.par,"\n")
  
  e0<-data[index,"Data1"]
  ea<-data[index,"Outside"]
  inv.cost<-data[index,"Investment1"]
  p0<-0.8
  p1<-0.5
  penalty<-data[index,"Penalty"]
  
#  par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,data[index,"old_attack_outcome"])
  
    old.payment.outcome<-ifelse(is.na(data[index,"old_payment_outcome"]),0.5,data[index,"old_payment_outcome"])

  # variation #2, averaged over t-1 and t-2
  
  # old.payment.outcome1<-ifelse(is.na(data[index,"old_payment_outcome"]),0.5,data[index,"old_payment_outcome"])
  # old.payment.outcome2<-ifelse(is.na(data[index,"old_t2_payment_outcome"]),0.5,data[index,"old_t2_payment_outcome"])
  
  # old.payment.outcome<-0.5*(old.payment.outcome1+old.payment.outcome2)
  
#  varation #1, using average in the past  
#  old.payment.outcome<-ifelse(is.na(data[index,"old_average_payment_outcome"]),0.5,data[index,"old_average_payment_outcome"])
  
  old.attack.outcome<-ifelse(is.na(data[index,"old_attack_outcome"]),0,data[index,"old_attack_outcome"])
#  old.ransom<-ifelse(is.na(data[index,"old_ransom_amount"]),0,data[index,"old_attack_outcome"])
  old.ransom<-ifelse(is.na(data[index,"old_ransom"]),0,data[index,"old_attack_outcome"])
  
  chat.invest<-data[index,"chat_invest"]
  chat.payment<-data[index,"chat_pay"]
  inv.social.norm<-ifelse(is.na(data[index,"old_i_1"]),0.5,0.5*(data[index,"old_i_1"]+data[index,"old_i_2"]))
  
  period<-data[index,"Period"]

#  variation #1,using average in the past  
#  inv.social.norm<-ifelse(is.na(data[index,"old_average_i_1"]),0.5,0.5*(data[index,"old_average_i_1"]+data[index,"old_average_i_2"]))

#  inv.social.norm.t1<-ifelse(is.na(data[index,"old_i_1"]),0.5,0.5*(data[index,"old_i_1"]+data[index,"old_i_2"]))
#  inv.social.norm.t2<-ifelse(is.na(data[index,"old_t2_i_1"]),0.5,0.5*(data[index,"old_t2_i_1"]+data[index,"old_t2_i_2"]))
  
#  inv.social.norm<-(1-final.par[10])*inv.social.norm.t1+final.par[10]*inv.social.norm.t2
#  inv.social.norm<-0.5*inv.social.norm.t1+0.5*inv.social.norm.t2
  
#  par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,old.attack.outcome,old.ransom,inv.social.norm)
#  par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,chat.invest,chat.payment,inv.social.norm)
  par.setting<-c(e0,ea,inv.cost,p0,p1,penalty,old.payment.outcome,0,0,inv.social.norm,period)
  
#  browser()
  
  gamma.i<-final.par[4]
  
  u.m1<-matrix(nrow=2,ncol=2)
  u.m2<-matrix(nrow=2,ncol=2)
  
  for(i in 0:1) for(j in 0:1) {
    u.m1[(i+1),(j+1)]<-u.inv(i,j,c(data[index,"old_i_1"],data[index,"old_i_2"]),par.setting = par.setting,par.behavior = final.par)
    u.m2[(i+1),(j+1)]<-u.inv(i,j,c(data[index,"old_i_2"],data[index,"old_i_1"]),par.setting = par.setting,par.behavior = final.par)
  }
  
  dimnames(u.m1)<-list(c("0","1"),c("0","1"))
  dimnames(u.m2)<-list(c("0","1"),c("0","1"))
  
#  if(sum(is.na(u.m1))>0) browser()
  
  qre.i<-find.qre(c(final.par[4],final.par[4]),list(u.m1,u.m2),tol=1e-8,max.iter=100,mode="m")

#if(final.par[4]<0.0001) browser()
  
  qre.i1<-qre.i[[1]]
  llk.i1<-(-log(qre.i1[qre.i1[,1]==data[index,"i_1"],2]))

  qre.i2<-qre.i[[2]]
  llk.i2<-(-log(qre.i2[qre.i2[,1]==data[index,"i_2"],2]))
  
  # browser()
  
  attack.decision<-data[index,"attack.decision"]
  
  qra<-qre.att(c(data[index,"i_1"],data[index,"i_2"]),par.setting=par.setting,par.behavior=final.par)
  llk.a<-(-log(qra[qra[,1]==data[index,"attack.decision"],2]))
  
  llk.r<-0
  r<-data[index,"ransom"]
  
  if(!is.na(r)) {
    
    if(r==0) r<-1
    if(attack.decision==1) inv<-data[index,"i_1"]
    if(attack.decision==2) inv<-data[index,"i_2"]
    
    qr<-qre.r(inv,par.setting=par.setting,par.behavior=final.par)
    llk.r<-(-log(qr[match(r,qr[,1]),2]))
#    browser()
  }
  
  llk.pay<-0
  if(!is.na(data[index,"p_1"])|!(is.na(data[index,"p_2"]))) {
    if(attack.decision==1) inv<-data[index,"i_1"]
    if(attack.decision==2) inv<-data[index,"i_2"]

    if(attack.decision==1) pay<-data[index,"p_1"]
    if(attack.decision==2) pay<-data[index,"p_2"]
    
    pp<-prob.pay(r,inv,par.setting=par.setting,par.behavior=final.par)    
    
    if(pp==1) llk.pay<-(1-pay)*100
    else llk.pay<-(-pay*log(pp)-(1-pay)*log(1-pp))
  }
  
  result<-llk.i1+llk.i2+llk.a+llk.r+llk.pay
  if(is.na(result)) browser()
  result
}


llk<-function(par,par.restrict=NULL,par.restrict.value=NULL,data=NULL) {
  
  # weird bug .. parallel gives different results!
  
result<-parLapply(cl=cluster,c(1:nrow(data)),llk.i,par=par,par.restrict=par.restrict,
                  par.restrict.value=par.restrict.value,data=data)
#  result<-lapply(c(1:nrow(data)),llk.i,par=par,par.restrict=par.restrict,
#                       par.restrict.value=par.restrict.value,data=data)
  #browser()
  #sum(result)
  do.call(sum,result)
}

fit.model<-function(start,data.list=NULL,lower=NULL,upper=NULL,scale=NULL,
                     par.restrict=NULL,par.restrict.value=NULL) {
  
   ans<-nlminb(start,scale=scale,lower=lower,upper=upper,llk,par.restrict=par.restrict,
               par.restrict.value=par.restrict.value,data=data.list,
               control=list(trace=2,iter.max=1000,eval.max=5000,rel.tol=1e-8))
#     a<-scale
#     ans<-optim(start,lower=lower,upper=upper,llk,par.restrict=par.restrict,
#                 par.restrict.value=par.restrict.value,data=data,method="L-BFGS-B",
#                 control=list(trace=6,maxit=1000))
    
    ans
}


